ZHAO Ang (赵昂)'s TA report for assignment03

SID: 12132606

Github: https://github.com/ZINSIST/ESE5023_Assignments_12132606

Responsible TA: HUANG Hao

Grade: 35

import numpy as np
import pandas as pd
import xarray as xr
from matplotlib import pyplot as plt
from matplotlib import cm
%matplotlib inline

1. Niño 3.4 index (7/10)

1.1

sst_Data = xr.open_dataset('NOAA_NCDC_ERSST_v3b_SST.nc')
sst_region = sst_Data.sel(lat=slice(-5,5), lon=slice(190,240))

group_data = sst_region.sst.groupby('time.month')
sst_anm = group_data - group_data.mean(dim='time')

sst_anm.values

array([[[-0.43157768, -0.41846275, -0.39795303, ..., -0.2116642 ,
         -0.23776245, -0.24401474],
        [-0.41259003, -0.4067192 , -0.3875141 , ..., -0.52064896,
         -0.5346451 , -0.51997185],
        [-0.40932274, -0.39743805, -0.36237717, ..., -0.6373882 ,
         -0.6171951 , -0.583725  ],
        [-0.4140854 , -0.37909317, -0.3215618 , ..., -0.43292618,
         -0.38404274, -0.3352623 ],
        [-0.5043678 , -0.43894005, -0.3710251 , ..., -0.17453575,
         -0.11044502, -0.06918144]],

       [[-0.5374584 , -0.52739716, -0.50823593, ..., -0.40254593,
         -0.44382668, -0.45287704],
        [-0.55093956, -0.539135  , -0.51673317, ..., -0.6660595 ,
         -0.7127285 , -0.710968  ],
        [-0.61242104, -0.5959244 , -0.5572338 , ..., -0.7235069 ,
         -0.7326374 , -0.73106194],
        [-0.6798363 , -0.6483364 , -0.5889931 , ..., -0.5397434 ,
         -0.50793266, -0.49977684],
        [-0.7830448 , -0.7286701 , -0.6683655 , ..., -0.33967972,
         -0.29167747, -0.27325058]],

       [[-0.4547863 , -0.43338776, -0.40325546, ..., -0.19753456,
         -0.23086166, -0.23381996],
        [-0.44275093, -0.40014458, -0.34834862, ..., -0.50234795,
         -0.53378105, -0.5277729 ],
        [-0.41475677, -0.34911728, -0.2700863 , ..., -0.6050377 ,
         -0.59763145, -0.59329414],
        [-0.4070778 , -0.33431816, -0.24388504, ..., -0.43722725,
         -0.40891075, -0.4024048 ],
        [-0.49510002, -0.42836   , -0.3584957 , ..., -0.23085403,
         -0.21210098, -0.20635033]],

       ...,

       [[-0.48802376, -0.55789375, -0.61621857, ..., -0.5144863 ,
         -0.446558  , -0.35310745],
        [-0.95103073, -1.0239639 , -1.0813885 , ..., -0.866415  ,
         -0.8075161 , -0.73119164],
        [-1.164526  , -1.2153225 , -1.2541466 , ..., -0.9210968 ,
         -0.87926865, -0.8485966 ],
        [-0.95791817, -1.0026321 , -1.020546  , ..., -0.5718498 ,
         -0.55397415, -0.55609703],
        [-0.48871994, -0.538929  , -0.56219673, ..., -0.34138107,
         -0.31702995, -0.2987461 ]],

       [[-0.55140495, -0.607254  , -0.6507721 , ..., -0.34615326,
         -0.2555828 , -0.13972664],
        [-0.989378  , -1.0497723 , -1.0954857 , ..., -0.86087227,
         -0.7690697 , -0.65498734],
        [-1.1887245 , -1.252285  , -1.3029232 , ..., -1.0460625 ,
         -0.9661274 , -0.8785801 ],
        [-1.002367  , -1.0756893 , -1.1325111 , ..., -0.7207298 ,
         -0.6597252 , -0.5900669 ],
        [-0.5770798 , -0.65514374, -0.72174263, ..., -0.4353485 ,
         -0.36265945, -0.28103828]],

       [[-0.3578701 , -0.41542053, -0.47110367, ..., -0.2400589 ,
         -0.1464405 , -0.03788376],
        [-0.7678585 , -0.83501625, -0.9024124 , ..., -0.727829  ,
         -0.61603355, -0.48027992],
        [-0.96187973, -1.0445309 , -1.1224213 , ..., -0.9327831 ,
         -0.81235695, -0.6655674 ],
        [-0.82112694, -0.9206734 , -1.0085506 , ..., -0.6531601 ,
         -0.5626869 , -0.4374504 ],
        [-0.4864292 , -0.5823746 , -0.6702862 , ..., -0.36221695,
         -0.30041504, -0.1987915 ]]], dtype=float32)

1.2

Your answer of question 1.2 is same as GONG Guoqing. 3 points were deducted.

sst_anm_rolling = sst_anm.rolling(time=3, center=True).mean()
line_anm = np.nanmean(sst_anm_rolling.values,axis=(1,2))

/var/folders/qy/ngsg7cz10ys7fn_5n46kknmr0000gn/T/ipykernel_68447/2672718069.py:2: RuntimeWarning: Mean of empty slice
  line_anm = np.nanmean(sst_anm_rolling.values,axis=(1,2))

time = pd.date_range(start='1960-01',periods=684,freq='m')
fig,ax = plt.subplots(1,1,figsize = [10,6],dpi=300)

ax.plot(time,line_anm,color='k')

#  xlabel, ylabel and title
ax.set_ylabel('Anomaly Degrees C', color='k', fontsize=15)
ax.set_xlabel('Year', color='k', fontsize=15)
ax.set_title("SST Anomaly in Niño 3.4 Region", fontsize=20)

ax.grid(linestyle='--',linewidth=0.3,alpha=0.5,color='k')

# put hlines into the figure
ax.hlines(y = 0.5,xmin=time[0],xmax=time[-1],color='r',linestyles='--',lw=0.5,label='EI Nino Threshold')
ax.hlines(y = -0.5,xmin=time[0],xmax=time[-1],color='b',linestyles='--',lw=0.5,label='La Nina Threshold')
ax.hlines(y = 0,xmin=time[0],xmax=time[-1],color='k',linestyles='solid',lw=1,label='3 mth running mean')

ax.set_ylim(-3,3)

ax.legend(loc='best',fontsize=12)

ax.fill_between(time,0,line_anm,where=(line_anm>0),color='r')
ax.fill_between(time,0,line_anm,where=(line_anm<0),color='b')

<matplotlib.collections.PolyCollection at 0x1124c2f40>

2. Earth’s energy budget (23/30)

2

ds_toa=xr.open_dataset('CERES_EBAF-TOA_200003-201701.nc')
da_lw=ds_toa['toa_lw_all_mon']
da_sw=ds_toa['toa_sw_all_mon']
da_sr=ds_toa['solar_mon']
da_nf=ds_toa['toa_net_all_mon']
time=ds_toa['time']

fig,axs=plt.subplots(2,2,sharex=True,sharey=True,figsize=(20,8))
da_lw.mean(dim='time').plot(ax=axs[0,0])
da_sw.mean(dim='time').plot(ax=axs[0,1])
da_sr.mean(dim='time').plot(ax=axs[1,0])
da_nf.mean(dim='time').plot(ax=axs[1,1])

axs[0,0].set_title('long wave',fontsize=15)
axs[0,1].set_title('short wave',fontsize=15)
axs[1,0].set_title('solar radiation',fontsize=15)
axs[1,1].set_title('net flux',fontsize=15)

Text(0.5, 1.0, 'net flux')

fig, ax = plt.subplots(1,1)
x = ds_toa.toa_sw_all_mon.lon
y = ds_toa.toa_sw_all_mon.lat
z = ds_toa.solar_mon.mean(dim='time') - (ds_toa.toa_lw_all_mon.mean(dim='time') + ds_toa.toa_sw_all_mon.mean(dim='time'))
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('ds_toa_solar - (lw + sw)')
cntr = ax.contourf(x, y, z, levels=30, cmap=cm.coolwarm)
fig.colorbar(cntr, ax=ax)
plt.show()

I think you did not verify (visually) that the addition is equivalent to the TOA net flux. I do it like follow.

fig, ax = plt.subplots(1,1)
x = ds_toa.toa_sw_all_mon.lon
y = ds_toa.toa_sw_all_mon.lat
z = (ds_toa.solar_mon.mean(dim='time') - ds_toa.toa_lw_all_mon.mean(dim='time') - ds_toa.toa_sw_all_mon.mean(dim='time')) - ds_toa.toa_net_all_mon.mean(dim='time')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('(ds_toa_solar - lw - sw) - net')
cntr = ax.contourf(x, y, z, levels=30, cmap='bwr')
fig.colorbar(cntr, ax=ax)
plt.show()

2.2

ds_toa=xr.open_dataset('CERES_EBAF-TOA_200003-201701.nc')
da_lw=ds_toa['toa_lw_all_mon']
da_sw=ds_toa['toa_sw_all_mon']
da_sr=ds_toa['solar_mon']
da_nf=ds_toa['toa_net_all_mon']

R=6371.4e3

lat_j=np.linspace(-89.5,89.5,180)
lat_j_rad=np.deg2rad(lat_j)

cos_j=np.cos(lat_j_rad)
Sj=2*np.pi**2*cos_j*R**2/(180*360)

Sij=np.repeat(Sj.reshape(-1,1),360,axis=1)
Stij=np.repeat(Sij.reshape(1,180,360),203,axis=0)

da_sr.values*=Stij
da_lw.values*=Stij
da_sw.values*=Stij

# the area of the whole surface of earth
area_glob=4*np.pi*R**2

sr=da_sr.mean(dim='time').values.sum()/area_glob
lw=da_lw.mean(dim='time').values.sum()/area_glob
sw=da_sw.mean(dim='time').values.sum()/area_glob

print('solar radiations (Wm-2):',sr.round(1))
print('long wave outgoing (Wm-2):',lw.round(1))
print('short wave outgoing (Wm-2):',sw.round(1))

solar radiations (Wm-2): 340.3
long wave outgoing (Wm-2): 240.3
short wave outgoing (Wm-2): 99.1

I find the 2.3 and 2.4 of your assignment03 is totally the same as CHEN Penghan, so I give you half points (5') for these two questions.

2.3

da_nf=ds_toa['toa_net_all_mon']
da_nf.values*=Stij/1e9

fig,axs=plt.subplots(1,2,figsize=(12,3),dpi=300)

da_nf.rename('').mean(dim='lon').transpose().plot(ax=axs[0],robust=True)
axs[0].set_title('Net radiation in each 1-degree latitude band $(10^9W)$')
da_nf.mean(dim={'lon','lat'}).plot(ax=axs[1])
axs[1].set_ylabel('net radiation $(10^9W)$')
axs[1].set_title('Net radiation in each 1-degree latitude band $(10^9W)$')

Text(0.5, 1.0, 'Net radiation in each 1-degree latitude band $(10^9W)$')

2.4

arrclda = ds_toa.cldarea_total_daynight_mon.mean(dim='time').values
high_cloud_area = (arrclda>=75)
low_cloud_area = (arrclda<=25)

hclw = ds_toa.toa_lw_all_mon.mean(dim='time')
hcsw = ds_toa.toa_sw_all_mon.mean(dim='time')
lclw = ds_toa.toa_lw_all_mon.mean(dim='time')
lcsw = ds_toa.toa_sw_all_mon.mean(dim='time')

# Create Figure and Subplots
fig,axes = plt.subplots(2,2,figsize=(16,8),dpi=500)
# Plot each axes
hclw.where(high_cloud_area).plot(ax=axes[0,0])
hcsw.where(high_cloud_area).plot(ax=axes[0,1])
lclw.where(low_cloud_area).plot(ax=axes[1,0])
lcsw.where(low_cloud_area).plot(ax=axes[1,1])
axes[0,0].set_title('high cloud area outgoing longwave',fontsize = 14)
axes[0,1].set_title('high cloud area outgoing shortwave',fontsize = 14)
axes[1,0].set_title('low cloud area outgoing longwave',fontsize = 14)
axes[1,1].set_title('low cloud area outgoing shortwave',fontsize = 14)
plt.tight_layout()

After our three TAs discussed about question 2.4 and 2.5. We think that's more proper to select radiation by cloud condition and then calculate mean values. That's because radiation values are affected by cloud condition at that moment, not by the mean conditions. So, modified your code as follow, please think more.

## By HUANG Hao
# Create Figure and Subplots
fig,axes = plt.subplots(2,2,figsize=(16,8),dpi=500)
# Plot each axes
ds_toa.toa_lw_all_mon.where(ds_toa.cldarea_total_daynight_mon >= 75).mean(dim='time').plot(ax=axes[0,0])
ds_toa.toa_sw_all_mon.where(ds_toa.cldarea_total_daynight_mon >= 75).mean(dim='time').plot(ax=axes[0,1])
ds_toa.toa_lw_all_mon.where(ds_toa.cldarea_total_daynight_mon <= 25).mean(dim='time').plot(ax=axes[1,0])
ds_toa.toa_sw_all_mon.where(ds_toa.cldarea_total_daynight_mon <= 25).mean(dim='time').plot(ax=axes[1,1])
axes[0,0].set_title('high cloud area outgoing longwave',fontsize = 14)
axes[0,1].set_title('high cloud area outgoing shortwave',fontsize = 14)
axes[1,0].set_title('low cloud area outgoing longwave',fontsize = 14)
axes[1,1].set_title('low cloud area outgoing shortwave',fontsize = 14)
plt.tight_layout()

2.5

You did not calcualted the global mean values and you also did not answer the question about what is the overall effect of clouds on shortwave and longwave radiation. 2 point was deducted.

print('high cloud long wave:',np.nanmean(hclw),'(W/m2)')
print('high cloud short wave:',np.nanmean(hcsw),'(W/m2)')
print('low cloud long wave:',np.nanmean(lclw),'(W/m2)')
print('low cloud short wave:',np.nanmean(lcsw),'(W/m2)')

high cloud long wave: 1891498600000.0 (W/m2)
high cloud short wave: 780467600000.0 (W/m2)
low cloud long wave: 1891498600000.0 (W/m2)
low cloud short wave: 780467600000.0 (W/m2)

3. Explore a netCDF dataset (5/10)

Your answer of question 3 is same as GONG Guoqing (龚国庆). 5 points were deducted.

3.1

ds = xr.open_dataset("air.sig995.2012.nc")

/Users/haohuang/miniconda3/envs/ese5023/lib/python3.9/site-packages/xarray/coding/times.py:123: SerializationWarning: Ambiguous reference date string: 1-1-1 00:00:0.0. The first value is assumed to be the year hence will be padded with zeros to remove the ambiguity (the padded reference date string is: 0001-1-1 00:00:0.0). To remove this message, remove the ambiguity by padding your reference date strings with zeros.
  warnings.warn(warning_msg, SerializationWarning)

group_data = ds.air.groupby('time.month')
air_anom = group_data - group_data.mean(dim='time')
line_air_anom = air_anom.mean(dim={'lat','lon'})

time = pd.date_range(start='2012-01-01',periods=366,freq='d')
fig,ax = plt.subplots(1,1,figsize = [10,6],dpi=300)

ax.plot(time,line_air_anom,color='k')

ax.set_ylabel('Anomaly Degrees C', color='k', fontsize=15)
ax.set_xlabel('Time', color='k', fontsize=15)
ax.set_title("2012 air temperature anomalies on the 0.995 sigma level", fontsize=20)


ax.grid(linestyle='--',linewidth=0.3,alpha=0.5,color='k')

ax.hlines(y = line_air_anom.max(),xmin=time[0],xmax=time[-1],color='khaki',linestyles='--',lw=0.5,label='Anomaly MAX positive value')
ax.hlines(y = line_air_anom.min(),xmin=time[0],xmax=time[-1],color='pink',linestyles='--',lw=0.5,label='Anomaly MAX negative value')
ax.hlines(y = 0,xmin=time[0],xmax=time[-1],color='k',linestyles='solid',lw=1,label='anomalies = 0')

ax.set_ylim(-2,2)
ax.legend(loc='best',fontsize=12)

ax.fill_between(time,0,line_air_anom,where=(line_air_anom>0),color='khaki')
ax.fill_between(time,0,line_air_anom,where=(line_air_anom<0),color='pink')

<matplotlib.collections.PolyCollection at 0x13d7749d0>

3.2

fig, axes = plt.subplots(2,3, figsize=(12,6), sharex=False, sharey=False, dpi=300)
da_air = ds.air
da_air_Dec = ds.air.sel(time=slice('2012-12-01', '2012-12-31'))
da_air_Jul = ds.air.sel(time=slice('2012-07-01','2012-07-31'))
da_air_shenzhen = ds.air.sel(lon='114',lat='22.5',method='nearest')
da_air_Ant = ds.air.sel(lat='-90',method='nearest')



da_air_Dec.mean('time').plot(ax=axes[0,0])
da_air_Jul.mean('time').plot(ax=axes[1,0],cmap='rainbow')
da_air_Dec.mean('lon').transpose().plot(ax=axes[0,1])
da_air_Jul.mean('lon').transpose().plot(ax=axes[1,1],cmap='rainbow')
da_air_shenzhen.plot(ax=axes[1,2],c='r')
da_air_Ant.mean('lon').plot(ax=axes[0,2])


axes[0,2].grid(linestyle='--', linewidth=0.5, alpha=0.5)
axes[1,2].grid(linestyle='--', linewidth=0.5, alpha=0.5)


axes[0,0].set_title('Mean temperature in Dec 2012 (K)',fontsize = 14)
axes[1,0].set_title('Mean temperature in Jul 2012 (K)',fontsize = 14)
axes[0,1].set_title('Temperature in Dec 2012 mean lon (K)',fontsize = 14)
axes[1,1].set_title('Temperature in Jul 2012 mean lon (K)',fontsize = 14)
axes[1,2].set_title('Mean temperature in ShenZhen 2012 (K)',fontsize = 14)
axes[0,2].set_title('Mean temperature in Antarctica 2012 (K)',fontsize = 14)

# better layout
plt.tight_layout()